#!/bin/bash

set -e
tools_path=/mnt/ilustre/app/medical/tools
data_path=/mnt/ilustre/app/medical/.data

source ${tools_path}/.var

r_path=${tools_path}/R


while getopts  "g:n:" opts
do
        case  $opts  in
        n)
			sample_name=$OPTARG
		;;
		g)
			ref_gene=$OPTARG
		;;
		\?)
			echo "Usage: $0 [-n sample_name] [-g ref_gene] <sample.bam> <ref.bam> <interval.bed>"
			echo -e "\t-n sample name (sap) "
			echo -e "\t-g ref gene file (refGene.txt)"
            exit 1
		;;
        esac
done
shift $(($OPTIND - 1))


if [ -z "$1" -o -z "$2" ]; then
	echo $0 "[-n sample_name] [-g ref_gene] <sample.bam> <ref.bam> <interval.bed>"
	exit 1
fi


script_path=$tools_path/script
quandico=${tools_path}/quandico-1.13/QUANDICO-v1.12/scripts/quandico
qgetcounts=${tools_path}/quandico-1.13/QUANDICO-v1.12/scripts/qgetcounts
qcluster=${tools_path}/quandico-1.13/QUANDICO-v1.12/scripts/qcluster
samtools=${tools_path}/samtools-1.2/samtools



if [ -z "$sample_name" ]; then
	sample_name=sap
fi

if [ -z "$ref_gene" ]; then
	ref_gene=${data_path}/refGene.txt
fi


# dat=`date +'%b_%d_%H_%M_%S_%Y'`

# dir_name=${dat}__${sample_name}
# mkdir $dir_name
# cd $dir_name

log=.log
if [ ! -e "$log" ]; then
	:> $log
fi

echo 2>>$log 1>&2
echo Sample name is: $sample_name 2>>$log 1>&2
echo $ref_gene 2>>$log 1>&2

echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo $1 2>>$log 1>&2
echo $2 2>>$log 1>&2
echo $3 2>>$log 1>&2

if [ ! -e ".Rprofile" ]; then
	cp $r_path/.Rprofile . 2>>$log
fi

###########
samp=$1
ref=$2
interval=$3

primerlen=22

###########


perl $qgetcounts \
-i $samp \
-a $interval \
-o 1.count.tsv \
--primerlen $primerlen \
--samtools $samtools \
--verbose 0 \
2>>$log

perl $qgetcounts \
-i $ref \
-a $interval \
-o 2.count.tsv \
--primerlen $primerlen \
--samtools $samtools \
--verbose 0 \
2>>$log


perl $qcluster \
-i 1.count.tsv \
--names $ref_gene \
> 1.cluster.tsv \
2>>$log


perl $qcluster \
-i 2.count.tsv \
--names $ref_gene \
> 2.cluster.tsv \
2>>$log



perl $quandico \
--no-cluster \
-A $interval \
-s data=1.cluster.tsv \
-r data=2.cluster.tsv \
-s x=2 -s y=0 \
-r x=2 -r y=0 \
--rexe $R \
--inject $r_path/quandico.r \
-d . \
-b 1 \
2>>$log


# echo 2>>$log 1>&2
# echo 2>>$log 1>&2
# echo quandico 2>>$log 1>&2
# perl $quandico \
# --counter $qgetcounts \
# --engine $qcluster \
# -A ${data_path}/intervals/1/brca.qiagen.NGHS-001X-Covered.b37.bed \
# -X samtools=$samtools \
# --rexe $R \
# -s map=$1 \
# -r map=$2 \
# -s x=2 -s y=0 \
# -r x=2 -r y=0 \
# --inject $r_path/quandico.r \
# --cp names=$ref_gene \
# -d . \
# -b 1 \
# 2>>$log

rm .Rprofile 2>>$log
